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ABSTRACT 


A two-step hybrid perturbation- Galerkin method to solve a variety of differential equations 
which involve a parameter is presented and discussed. The method consists of: (1) the use 
of a perturbation method to determine the asymptotic expansion of the solution about one 
or more values of the parameter; and (2) the use of some of the perturbation coefficient 
functions as trial functions in the classical Bubnov-Galerkin method. This hybrid method 
has the potential of overcoming some of the drawbacks of the perturbation method and the 
Bubnov-Galerkin method when they are applied by themselves, while combining some of 
the good features of both. The proposed method is illustrated first with a simple linear two- 
point boundary value problem and is then applied to a nonlinear two-point boundary value 
problem in lubrication theory. The results obtained from the hybrid method are compared 
with approximate solutions obtained by purely numerical methods. Some general features 
of the method, as well as some special tips for its implementation, are discussed. A survey 
of some ^current jresearch application areas is presented and its degree of applicability to 
broader problem areas is discussed. 


1 This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS 1-18605 while the first author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 


The method we wish to describe is a two-step hybrid analysis technique, which was ap- 
parently first studied by Ahmed K. Noor and collaborators in conjunction with the finite 
element analysis of geometrically nonlinear problems in structural mechanics (see Geer 
and Andersen (1989a) for several references). It is based upon the successive use of per- 
turbation expansion methods and the classical Bubnov-Galerkin approximation technique. 
In the perturbation method, an approximation to the solution of a particular problem 
involving a parameter is developed in terms of a series of unknown functions with preas- 
signed coefficients, i.e. gauge functions. The unknown functions are determined by solving 
a recursive set of differential equations which are, in general, simpler than the original 
governing differential equation. By contrast, in the Bubnov-Galerkin technique one seeks 
an approximate solution to the problem in the form of a linear combination of specified 
(known) coordinate functions with unknown coefficients. The coefficients are determined 
by demanding that the residual formed by substituting the trial solution into the governing 
differential equation is orthogonal to each of the coordinate functions. 

The Galerkin method (see Galerkin (1915)) has, of course, been known and used for 
a long time. But a principle problem associated with i+s successful application lies in 
the choice of appropriate basis functions. In a series of papers Noor and his collaborators 
(Noor et al (1979-1988)) have shown for a variety of structural mechanics problems that the 
terms in a Taylor series expansion of the solution of a parameterized system of discretized 
equations can be particularly effective as Galerkin trial functions (or basis vectors). In 
addition, it has been repeatedly demonstrated that the “reduced-basis” solutions can be 
useful for significantly larger values of the expansion parameter than are the Taylor series 
solutions on which they are based. Noor and Balch (1984) and Noor, Balch and Shibut 
(1984) have also applied the same general principles but without discretization to some 
thermal analysis problems. A treatment of the reduced basis method from a mathematical 
point of view is given by Fink and Rheinbolt (1983). 

Some general observations about the technique are the following. First, in many per- 
turbation problems, much effort has to be expended to compute each additional term in 
a perturbation expansion. Through the use of the proposed hybrid method, the known 
perturbation terms can be exploited more fully. Secondly, another way of viewing the 
technique is to recognize that in many perturbation expansions the functional form of the 
higher-order terms can be well approximated by a linear combination of the lower order 
terms. Thus, much of the effect of the higher-order terms may be included by applying 
the reduced basis technique to the lower order terms. Finally, our investigations indicate 
that, while the use of a Taylor series expansion is frequently limited by a finite radius of 
convergence, the proposed hybrid method can sometimes yield good results well outside 
the radius of convergence. This is true even if that radius is zero (see Andersen and Geer 
(1988)). 

It is our belief that the junction of perturbation and Galerkin techniques can be useful in 
a wide variety of application areas and in these applications the hybrid technique will give 
better approximations than the perturbation method alone. In this paper we will present 
applications of the technique independent of finite element or finite difference methods and 
in an area apart from structural mechanics. In particular, in the next section, we shall 
describe the method in more detail and then illustrate it with a simple linear two-point 
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boundary value problem in section 3. In section 4 we apply the method to a nonlinear 
problem in lubrication theory and use the results to determine the load carrying capacity 
of a slider bearing. In the final section, we make some general observations and comments 
on the method, as well as survey some of the current research application areas. 

2. DESCRIPTION OF THE METHOD 

The method we wish to describe is a two-step hybrid analysis technique. It is based 
upon the successive use of a straightforward perturbation expansion method and the classi- 
cal Bubnov- Galerkin approximation technique, as outlined in the introduction. While each 
of these methods is useful and has been successful in providing approximate solutions to a 
wide variety of nonlinear (and otherwise difficult) problems, each of these techniques has 
certain drawbacks. The perturbation method has at least two major drawbacks. First, as 
the number of terms in the perturbation expansion increases, the mathematical complexity 
of the equations which determine the unknown functions increases rapidly. Thus, in most 
practical applications, the perturbation series is limited to only a few terms. A second 
drawback to the perturbation method is the requirement of restricting the perturbation 
parameter to small values in order to obtain solutions of acceptable accuracy. (These 
drawbacks of the perturbation method have been recognized and several modifications or 
extensions have been proposed, see e.g. Van Dyke (1974) and Andersen and Geer (1982).) 
The main shortcoming of the Bubnov-Galerkin method is the difficulty, from a practical 
point of view, of selecting good coordinate functions. 

To illustrate the general ideas of the hybrid (or “reduced basis”) method, suppose we 
are seeking (an approximation to) the solution u(x) to the problem 

C(u,e) = 0, a: in some domain (or interval) V, (1) 

where C is some differential operator and e is a parameter. Without loss of generality, we 
can assume that the boundary conditions are homogeneous in u. We first generate the 
coordinate functions in a perturbation expansion of u about one or more specific values of 
the parameter e, say about e = e p , p = 1, 2, . . . , P, in the form 

n,-l 

U = £ K a *( f ) + °( a n r ( € ))> ( 2 ) 

k=0 

where {a*(e)} is an appropriate asymptotic sequence of gauge functions (e.g., a\ = e k 
for Taylor series expansions) and each u\ can be determined completely by a standard 
perturbation method (e.g. a composite expansion of inner and outer expansions). A 
subset of all of the perturbation functions u p k are now chosen as coordinate functions for 
the Bubnov-Galerkin technique and an approximation u for u is sought in the form 

N 

U = Y1 6j Uj , (3) 

;=1 

where the (unknown) parameters {6, } are functions of e and represent the amplitudes of 
the coordinate functions Uj . Here each Uj is one of the perturbation coordinate functions 
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u r k . To determine the unknown amplitudes { <5 ; } , we apply the Bubnov-Galerkin technique 
to the governing equation (1). Thus, we substitute (3) into (1) and demand that the 
residual be orthogonal to the N coordinate functions over the domain T>, i.e. 


f N 

/ C(%2 6, v,j,e)u k doc = 0, k = 1 , 2 ,..., AT. (4) 

* ;=i 

Equations (4) represent a set of N equations for the N unknown amplitudes {<5 ; }. While 
equations (4) must, in general, be solved numerically, solving them is simpler than nu- 
merically solving (1). In particular, for a fixed value of e, the solution to (4) is a point in 
N- dimensional space, where N is reasonably small, while the solution of (1) is a continuous 
function. 

We should note that this particular choice of coordinate functions overcomes the main 
drawback of the Bubnov-Galerkin method. By the way they are constructed, the pertur- 
bation coordinate functions are (under certain assumptions) elements of a set of functions 
which span the space of solutions in a neighborhood of their point of generation. Thus, 
they should fully characterize the solution u in that neighborhood. Also, in many applica- 
tions, the functions u k are determined by solving a set of linear equations, even though the 
original operator C may be nonlinear. The first property is necessary for the convergence 
of the Bubnov-Galerkin method, while the second property enhances the effectiveness of 
the proposed hybrid method for solving nonlinear problems. 

Another important property of the proposed method is that the coordinate functions, 
i.e. the perturbation functions, do not need to come from a regular perturbation expansion. 
In fact, all that is needed is a formal asymptotic expansion of the solution to (1) for values 
of e close to e p in the form of (2), where the (a*(e)} are a set of appropriate gauge 
functions. Here we are using the term “gauge functions” (Van Dyke (1975)) in the sense 
that the expressions {a£ } form an asymptotic sequence as e — * e p , i.e. for each fixed value 
of k the ratio /<** — » 0 as e — ► e p . Thus these expressions may involve terms like 
loge (Geer and Andersen (1989a)) or fractional powers of e (Geer and Andersen (1989b)). 
Hence, it follows that the proposed method may be applied to singular, as well as regular, 
perturbation problems. (See Geer and Andersen (1989a). In the following section, we shall 
illustrate our method using the regular perturbation expansion about e = 0 for a simple 
two-point boundary value problem, while in section 4 we shall apply the method using 
both a regular and a singular perturbation expansion of the solution. 


3. A SIMPLE EXAMPLE 

To illustrate the ideas just discussed, we consider the following simple two-point bound- 
ary value problem: 


u" + fu = e, 0 < x < 1, 
with u(0) = 0 and u(l) = 0. 
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The exact solution to this problem for e > 0 is 


u = 1 — cos(wx) 


1 — cos(w)' 
sin(w) 

u> = \fk. 


sin(wx), 


( 6 ) 


This function, which is graphed in Fig. 1 , has (real) singularities at e = 7r 2 , (3 tt) 2 , (5tt )^, . . .. 
A regular perturbation expansion of u yields a series of the form 

N 

u = 53 u A x ) ^ + °( cJV+1 )) (?) 

;=1 

where each Uj is a polynomial in x of degree 2 j which vanishes at x = 0 and x = 1. 
Following the ideas discussed above, we look for an approximate solution u in the form 

N 

u = J2 u i > ( 8 ) 

1 

where the amplitudes {6j } are to be determined for each value of e of interest. To determine 
them, we substitute (8) into (5) and apply the Bubnov-Galerkin criterion (4) to obtain the 
set of equations 


N 


a k,j(e) = M e )> & = 1,2,..., AT, 

j= 1 


( 9 ) 


where 


a* 


j(e) = u k (x)[u"( x ) + euj(x)]dx, 

b k {e) = e u k (x)dx. 

J r\ 


( 10 ) 


In this case the system of equations for the } is linear since the original problem is 
linear. For example, setting N = 2 in (9), we find 


6 1 = 


fa = 


1008e - 112e 2 
1008 - 112e + e 2 ’ 
1008e 2 

1008 - 112e + f 2 


( 11 ) 


For small values of e, we see that these expressions reduce to e and e 2 , respectively, as they 
should. Also, they remain finite as e becomes large, which suggests that our approximate 
hybrid solution might be valid for larger values of e than the perturbation solution on which 
it is based. For N = 4, we find by using the symbolic computation system Mathematica 
(Wolfram (1988)) that 


5 



( 12 ) 


<5i = —1232e(— 159120 + 18720c - 273c 2 + c 3 )/D, 

S 2 = 48048c 2 (4080 - 480c + 7 e 2 )/D, 

63 = — 11531520c 3 (— 17 + 2 e)/D, 

64 = 196035840 e 4 /D, 

D = 196035840 - 23063040c + 336336c 2 - 1232c 3 + c 4 . 

For small values of e we see that reduces to c J + 0(e 5 ). 

In Figs. 2 through 5 we compare two- and four- term hybrid solutions (8) with two- 
and four-term perturbation solutions (7) and with the exact solution (6). Fig. 2 shows 
comparisons for c = 2, 4 and 6. This figure indicates that the two-term perturbation 
solution ( N = 2) does well for e up to about c = 2, but that by c = 6 even the four-term 
e = 2, but that by e = 6 even the four-term perturbation solution (N=4) has serious errors. 
In this figure the two-term hybrid solutions agree very well with the exact solutions. 

Fig. 3 shows similar comparisons for c near the radius of convergence, 7r 2 , of the pertur- 
bation expansion. The two- and four-term perturbation solutions are not sensitive to the 
small variations in c and are poor approximations. The two-term hybrid solution is still 
very accurate. However, for values of c much closer to tt 2 than the values used in Fig. 3, the 
relative error of the hybrid solution becomes large because the two-term hybrid solution 
predicts the singularity to be at e = (3.141616) 2 rather than at tt 2 . Higher-order hybrid 
solutions become more accurate in this respect. The four-term hybrid solution gives the 
first singularity accurate to 13 digits. 

Fig. 4 compares the two- and four-term hybrid solutions with the exact solution at 
c = (27r) 2 . The two-term hybrid solution is a good approximation, but the four-term 
hybrid solution is much better. For values of c this large the perturbation solutions are 
useless. Finally, in Fig. 5 we compare for e = 1.02 (37t) 2 the two- and four-term hybrid 
solutions with the exact solution. The two-term hybrid solution is no longer useful, but 
the four-term hybrid solution is still satisfactory. In general, the larger e becomes the more 
terms are needed for the hybrid solution to give a useful approximation. 

Some of the reasons why the hybrid method provides a good approximation in this case, 
as well as insights this example provides for the validity of the method, will be discussed in 
section 5. Applications of the hybrid method to some general classes of two-point boundary 
value problems have been discussed in Geer and Andersen (1989b). 

4. APPLICATION - A SLIDER BEARING PROBLEM 

We now illustrate the hybrid method for a nonlinear example. We consider the problem 
of determining the excess pressure u(x) under a slider bearing, discussed for large bearing 
numbers by DiPrima (1969). Here u(x) satisfies 

4^ 3(I) 11 + “ <X)I 5^ = Tx {h{x)[1 + “ <X)IK (13) 

0 < x < 1, with u(0) = 0 = it(l). 

Here e, the inverse of the bearing number, is defined by e = p a d 2 /(6p B U), where p a is 
the ambient pressure in the fluid, d is the clearance of the trailing edge of the bearing, p is 
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the dynamic viscosity of the fluid, B is the length of the bearing in the flow direction, and 
U is the velocity of the bearing. The function h(x) is the prescribed film thickness under 
the bearing normalized such that h( 1) = 1. 

For small values of e (large bearing numbers), DiPrima constructed an asymptotic 
expansion of the form (2) with e\ — 0 and aj = e j , j = 0 , 1 ,..., i.e. 

u — Uq + u\ e + u\e 2 + . . . (14) 

with 


“0 = rTTT + «— ) - MO) 


h(x) 


-(1 - x) {Q(i) - h( 0) + 1}, 


(15) 


, I h(0) — 1 — Q \ „ 1 — x 

h{0) log ! h(o) - 1 \ + Q ~ — 7~' 


For large values of e, u has a regular perturbation expansion of the form (2) with 


aj = *, i.e. 


u = Uq€ 1 + u\e 2 + u|e 3 + (16) 

Substituting (16) into (13) and equating coefficients of like powers of e -1 on each side 
of the resulting equation, we find that each u 2 satisfies a simple second-order differential 
equation with homogeneous boundary conditions. In particular, we find that uq satisfies 


whose solution is 


4~{h 3 ( x ) = J- h ( x )> u o(°) = 0 = u o(!)> 

CIX CIX UJL 


2 = / h 2 (s)ds + c f h 3 (s)ds, 
J 0 u 

o = —J h~ 2 (s)ds j J h~ 3 (s)ds. 


( 17 ) 


(18) 


0 "0 

We now select one term from the small-e expansion (14) of u and one term from the 
large-e expansion (16) and construct a hybrid solution u in the form (3) with N — 2, i.e. 


u = 6iul(x,e) + 62 Uq( x ). (19) 

The quadratically nonlinear equations (4) for the amplitudes {fy } become 

• f () 1 { e £[ /l3 [ 1 + “ 1 ^ ~ + w ]} }«(}(*> < 0 ^ = °> 

( 20 ) 

+ «]£] - + u}})ul{x)dx = 0. 

These equations can be solved easily (and efficiently) by using Newton’s method, starting 
at small values of e and then proceeding to larger values of e. 
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To illustrate some results obtained by this approach we consider first a wedge-shaped 
slider bearing with h(x ) = 2 — x. [For each of the figures in this section hybrid results 
are indicated by solid lines, small-e perturbation expansions by lines with short dashes, 
large-e perturbation expansions by lines with long dashes, and the numerical results (which 
were obtained by a shooting method) with circles.] In Fig. 6, we compare for e = 0.7 the 
hybrid solution for N = 2 with the one-term small-parameter solution (14), the one-term 
large-parameter solution (16) and the solution obtained by a numerical shooting method. 
The value 0.7 was chosen for e because for the given function h(x) of all positive values 
of e this represents nearly the poorest agreement, an L 2 relative error norm of about 1%, 
between the hybrid solution and the numerical solution. A comparison of the hybrid and 
numerical solutions for four other values of e is shown in Fig. 7. In Fig. 8, we have plotted 
the amplitudes 61 and 62 (the coefficients of the small- and large-e coordinate functions, 
respectively) as functions of e. As the figure illustrates, 6\ approaches 1 and 62 approaches 
0 as e approaches zero, while 62 dominates <5i as e approaches infinity. In Figs. 9 and 10 
we have made comparisons similar to that shown in Fig. 6, this time for film thicknesses 
corresponding to a parabolic shape with h(x) = 1 + 2(1 — x) 2 , a case discussed by DiPrima, 
and to a cubic shape with h(x) = 1 + (1 — x) 2 (1 + 2 x). Again, the agreement between the 
hybrid and numerical solutions is very good, even for the “nearly worst case” values of e 
depicted in the figures. 

As an application of these results, we compute the nondimensionalized load carrying 
capacity W of the bearing, defined by IF = W /(p a B), where W is the load carrying 
capacity per unit width. W is given by 



Inserting our two- term hybrid approximation for u into (21), our approximation W to 
W becomes 


W = 6 1 W\ + 62 W 2 , 


( 22 ) 


where 


W 1 



= |(^(0) - 3) 


2*5(~)[ 1 + 2e + e( 2(“)l 



w 2 = 



x) dx 


= — f x[h 2 (x) + ch 3 (x)] dx, 


(23) 


(24) 
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and Uo(z,e) and ujj(x) are defined by (15) and (18), respectively. In obtaining (23) we 
have used some of the ideas of DiPrima to evaluate explicitly some of the integrals which 
appear in this expression. 

To compute W from formulas (22) - (24) for a specific bearing configuration, the 
prescription is: 

1) specify the film thickness profile h(x) and the bearing number c; 

2) calculate the quantity Q{ 1/e) from (15b) with x — 0; 

3) evaluate the integrals involving h(x) in (23) and (24) and use these results together 
with Q( 1/e) to find the values of W\ and W2; 

4) calculate the values of 8 \ and 82 from equations (20); and 

5) insert the values of 81, 82, W\, and W2 into (22) to calculate W . 

For the special case of a wedge slider bearing, where h(x) = 1 + «(1 — x) with a > 0, 
these equations become 


W\ = \{a - 2) + 1 + 0 log(l + o) 

2 a 

-lg(i){i + 2 £ + e <3(i)}, 

W 2 = -ilog(l + a)-— 2 

a 2 v ' a(2 + a) 


(25) 


Here Q( 1/e) = Q is determined from equation (15b) with x — 0, which thus becomes 


(1 + a) log [ — — | + Q = — — • (26) 

1 a 1 e 

In Fig. 11 we have plotted the load carrying capacity, W, for the wedge-shaped bearing 
with h(x) = 2 — x (i.e., a = 1) as e varies between 0 and 0.6. Shown are the two- 
term hybrid approximation, W ; the one-term small-e approximation, W\ (see Diprima 
(1969)); the one-term large-e approximation, which is W2/e; and a few points calculated 
numerically. As expected, the hybrid and the small-e perturbation calculations agree near 
e = 0. Also the hybrid and the large-e perturbation calculations agree for very large e. 


the large-e perturbation calculations agree for large e. While there is an intermediate 
range of values of e for which the small- e and large-e perturbation results agree very well, 
neither perturbation result has much accuracy there. On the other hand, it appears that 
the hybrid and numerical results are in very close agreement with each other throughout 
the entire range 0 < e < 00 . 


5. DISCUSSION 

We believe that the general technique we have outlined here provides a very useful 
method for analytically analyzing the behavior of solutions to differential equations over 
much wider ranges of parameter values than can be obtained by the perturbation or Galer- 
kin methods alone. The key seems to be that the perturbation functions are particularly 
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effective as “trial” functions for the Galerkin method, since they characterize the solution 
well for certain “extreme” values of the parameter. 

Some observations about the hybrid method, as we have described it, are the following. 
When higher-order terms in a perturbation expansion of the solution can be well approx- 
imated by linear combinations of lower order terms, the hybrid method appears to be an 
effective way to include some of the effects of these higher-order terms without actually 
computing them. In a practical problem, this property of the higher-order terms may 
manifest itself in the following way. A sequence of perturbation approximations based on 
lower order terms may appear to have approximately the right “shape” for the solution, 
but not the correct amplitude. This property is illustrated in our simple example in section 
3 for small values of e by comparing the perturbation approximations based on two and 
four terms. In this case, the hybrid approximation based on just the first two perturbation 
coefficient functions provides the necessary corrections for the “amplitudes” of these terms 
and results in a greatly improved approximate solution. More formally, this property can 
be expressed by the degree of linear dependence among the perturbation coordinate func- 
tions. If the perturbation functions are highly dependent, then only a few of them are 
needed in the hybrid solution to provide a reasonable approximation to the actual solu- 
tion. An extreme case of this was discussed by Geer and Andersen (1989a) in connection 
with a problem in slender body theory. In this case all of the perturbation coordinate 
functions are multiples of the first coordinate function and hence all of the functions are 
dependent. While the partial sums in the perturbation solution converge (slowly) to the 
exact solution, the hybrid solution based on only one perturbation coordinate function 
reproduces the exact solution to the problem. 

In a particular problem, the degree of dependence of a set of perturbation coordinate 
functions on the domain of interest may have to be obtained empirically. For example, 
the perturbation coordiante functions associated with our simple example in section 3 are 
polynomials of different orders and hence are definitely linearly independent. However, on 
the interval [0,1] they have a very high degree of linear dependency, which can be illustrated 
by normalizing the first several coordinate functions and plotting them on the same graph. 
(We have plotted the first fifteen (normalized) polynomials associated with our simple 
example and have found that they all have essentially the same parabolic shape, with 
only two or three curves being distinguishable at usual plotting resolution.) As a practical 
matter, a set of perturbation coordinate functions that are highly dependent may cause 
some numerical sensitivities when solving equations (4) for the amplitudes 6 , . Intuitively, 
this follows from the form of equations (4), which differ from each other only by the choice 
of the “test” functions u k which multiply the residual in the integrand. In this case, we 
have found it helpful to apply an orthonormalization process (such as the Gram-Schmidt 
process) to the the set {u*} and use the resulting set of orthonormal functions as the basis 
for the hybrid method. These steps will lead to the same hybrid approximation as would 
be obtained without the orthonormalization process, but the equations (4) will be better 
conditioned. 

It may be observed that from the Galerkin point of view the perturbation coordinate 
functions (7) generated for the simple two-point boundary- value problem are not very 
exciting. For example, the first two such functions are 
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ui = — x (1 — x)/2, 

U 2 = — x (1 — x) ( x 2 — x — l)/24. 


( 27 ) 


The factors x and 1 — x are needed to satisfy the boundary conditions; the problem has 
built-in symmetry about x = 1/2; and the overall numerical factors are immaterial. It 
turns out that the basis for the space spanned by the first N functions Uj is simply the 
sequence {:r(l — x)(x — l/2) 2 ^ -1 ^} with j = 1,2 ,...,jV, and this set of trial functions 
could easily have been conjectured for use in a stand-alone Galerkin solution. Nevertheless, 
it still may be said that the perturbation coordinate functions are good trial functions, 
that a very significant improvement over the perturbation result has been demonstrated, 
and that this improvement can be expected even for problems for which no “obvious” set 
of trial functions can be conjectured. 

In the simple example of section 3, we illustrated the fact that the hybrid method may 
produce meaningful approximations to the actual solution even well beyond the radius 
of convergence of the perturbation series. This is also illustrated by the simple example 
discussed by Geer and Andersen (1989a), as well as by the eigenvalue problem discussed by 
Andersen and Geer (1988). In the eigenvalue problem, the perturbation series about e = 0 
has a zero radius of convergence and hence is essentially useless for obtaining quantitative 
information for nonzero values of e. Hybrid solutions based on just a few terms in this 
series, however, provide a sequence of approximations which apparently converge to the 
correct solution for any e as the number of terms used increases. 

In addition, we have also found the hybrid method to be useful when dealing with a 
perturbation series with poor convergence properties within its radius of convergence. This 
was illustrated by Geer and Andersen (1989a) for a series which arises in the computation 
of the electrostatic potential field about a slender body of revolution. Here the perturbation 
solution is a series which contains an infinite number of terms, each one of which is an 
infinite series involving inverse powers of log e. These series appear to converge for small 
values of e of interest, but the rate of convergence is very slow. We were able to obtain a 
very meaningful approximation to the solution for these values of e by using just a few of 
the perturbation coordinate functions in the hybrid solution. 

In general, the amplitudes 6j which appear in the hybrid solution need to be determined 
numerically. The existence of relatively simple analytical expressions for the 6j , such as 
obtained in section 3 are definitely atypical. However, the availability of such expressions 
for a few simple examples provides us with useful insight into the method. 

As was pointed out in section 3, the <5 ; - for our simple example have real singularities 
(as functions of e) which approximate some of the singularities of the actual solution. (It 
is often these singularities which limit the convergence (and hence the range of usefulness) 
of perturbation solutions.) This property of the 6, seems to be a key reason why the 
hybrid solution can approximate functions (i.e., solutions to particular problems) which 
have various types of singularities. The solution of the simple example discussed by Geer 
and Andersen (1989a) has singularities corresponding to purely imaginary values of e. 
In this case, the corresponding 6j also have purely imaginary singularities, which again 
approximate some of the singularities of the solution. 

For the example discussed in section 4, it was relatively straightforward to obtain two 
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perturbation approximations to the solution - one valid formally for small values of e and 
the other valid formally for large values of e. (This property is also illustrated for several 
other classes of two-point boundary value problems for ordinary differential equations in 
Geer and Andersen (1989b).) In any particular problem, only one of these solutions may 
be physically relevant. However, if they can both be constructed, it is our recommendation 
that this be done (although, in practice, it may be much easier to compute higher-order 
terms in one expansion than in the other). The reason for this is that, in general, it has 
been our experience that it is better to use at least one term from each of these expansions 
in the hybrid approximation, than to use several terms from just one of the expansions 
(see especially the eigenvalue example discussed by Geer and Andersen (1989b).) 

It is also worth noting that only relatively simple mathematical operations are required 
to implement the hybrid method. In particular, once the perturbation coordinate functions 
have been determined, the Galerkin conditions (4) can be obtained by using a simple 
quadrature formula, such as Simpson’s rule. Once these equations are determined, they 
can be solved (if necessary) by Newton’s method. In practical terms, the Jacobian matrix 
required for this method can be well approximated by using a simple central difference 
formula and two more numerical quadratures for each component. 

Finally, we comment that if the differential equation to be solved has no parameter, 
then one should be invented. The given equation should be rewritten as an equation with 
a parameter. For one value of the parameter the new equation should coincide with the 
given equation. For another value of the parameter, the new equation should reduce to an 
equation which somehow is simpler than the original equation, perhaps by virtue of being 
linear instead of nonlinear, or by having higher symmetry. Then the asymptotic expansion 
is performed about this simpler equation. 

6. EXTENSIONS 

Numerous variations on the hybrid technique are possible. Some, such as expansions 
in two or more parameters simultaneously through the use of higher-order mixed partial 
derivatives have already found application. Other variations of the hybrid method need 
further investigation. For instance, the Bubnov-Galerkin step uses the perturbation coor- 
dinate functions as both trial and test functions. Perhaps, for some equations or systems 
of equations, a different weighted residual technique should be used in which the set of test 
functions should be different from the set of trial functions. It is known, however, that for 
a wide variety of problems the Bubnov-Galerkin step is simple in concept, effective, and 
easy to implement. 

The hybrid technique is effective for numerical problems as well as for problems which 
may be analyzed in terms of analytic functions. Many problems are far too complex 
because of their geometry, discontinuities or other reasons for effective analysis in terms 
of simple analytic functions. For these problems finite element, finite difference, Rayleigh- 
Ritz or other discretization techniques are commonly employed. However, the generation 
of series expansions and the formulation of the Galerkin step for discretized problems have 
direct numerical counterparts to the analytic steps employed in the study presented here. 
For discretized problems basis vectors in a large but finite-dimensional space take the place 
of the test and trial functions. The trick is to pick a small number of appropriate basis 
vectors to use as input to the Galerkin step, hence the term “reduced-basis method”. 
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We are currently investigating several possible application areas for the hybrid method 
and many of our preliminary results have been encouraging. In particular, we are investi- 
gating the application of the method to several eigenvalue problems, both for ordinary and 
partial differential equations. In addition to the quantum mechanical eigenvalue problem 
discussed in Andersen and Geer (1988) and Geer and Andersen (1989b) we are investi- 
gating some eigenvalue problems for partial differential equations where the boundary of 
the spacial domain for the problem is a perturbation of a “simple” domain for which the 
eigenvalue problem can be solved in closed form. 

We also feel that the method is potentially very useful for solving a variety of exte- 
rior boundary value problems for partial differential equations. In particular, we have 
already applied the method to some model problems in slender body theory (Geer and 
Andersen (1989a)) with promising results, and are currently working on combining the 
method with some simple homotopy ideas to further increase the accuracy of slender body 
approximations . 

The field of free and forced nonlinear vibrations is another application area we are inves- 
tigating. In this area we have found the method to be particularly effective for computing 
resonant frequencies of certain nonlinear oscillators (Geer and Andersen (1989b)), and are 
also currently using the method in conjunction with the multiple time scales method to 
compute frequency response curves for certain forced nonlinear oscillation problems. 
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Figure 1. Exact solution for simple diff 
of x and uj — v/(e) 




Figure 2. Comparison for the simple equation (5) of exact solutions (6), perturbation 
solutions (7) with N — 2 and N = 4, and hybrid solutions (8) with N = 2. 
The parameter e has values 2, 4, and 6. 
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Figure 3. Comparison for the simple equation (5) of exact solutions (6), perturbation 
solutions (7) with N = 2 and N = 4, and hybrid solutions (8) with N = 2. 
The parameter e has values (0.97 t) 2 , (0.977t) 2 , and (1.037 t) 2 , and (l.l7r) 2 . The 
perturbation solutions are indistinguishable from one another. 
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Figure 4. Comparison for the simple equation (5) of exact solutions (6) and hybrid so- 
lutions (8) with N = 2 and N = 4. The parameter e has the value (2.7r) 2 . 
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Figure 6. Comparison of solutions of the wedge-shaped slider bearing equation (13). 

One-term small-e perturbation solution (14), one-term large-e perturbation 
solution (16), two-term hybrid solution (19), and shooting-method (numerical) 
solution. 
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Figure 7. Comparison of pressure distributions for a wedge-shaped bearing with e 
0.1, 0.2, 0.3 and 1.0 as computed by hybrid and numerical methods. 
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